program heterogeneity_surfaces_ordered

    implicit none
    
    interface
     subroutine listing(SPACE, max)
     integer :: SPACE, max
     end subroutine
    end interface
    
    interface
     subroutine clinoform_filename(clinoform_filename_P, P, max, B, filepath)
     character :: clinoform_filename_P*1000,B*5, filepath*1000
     integer :: P, max
     end subroutine
    end interface

    interface
     subroutine clinoform_surfaces(TEMP,CLINOFORM_P,passage,clinoform_filename_P, n, m, P, nonvalue,xmin,xmax,ymin,ymax)
     character :: clinoform_filename_P*1000
     integer :: n, m, P
     real :: nonvalue,xmin,xmax,ymin,ymax
     real, allocatable :: CLINOFORM_P(:,:,:)
     real, allocatable :: TEMP(:,:,:)
     real, allocatable :: passage(:,:)
     end subroutine
    end interface
    
    interface
     subroutine more_nodes (BIGbis, nbis, mbis, n_regridbis, m_regridbis, NODE_CHANGEbis,Nonvalue, numsurf)
     integer :: n_bis,m_bis, n_regridbis, m_regridbis, numsurf
     real :: Nonvalue
     real, allocatable :: BIGbis(:,:)
     real, allocatable :: NODE_CHANGEbis(:,:)
     end subroutine
    end interface
    
    interface
     subroutine listing_surfaces(max,A,P)
     character :: A*5
     integer :: max, P
     end subroutine
    end interface
    
    interface
     subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     character :: temp_name*1000
     integer :: n,m, k
     real :: xmin,xmax,ymin,ymax
     real, allocatable :: TEMP(:,:,:)
     end subroutine
    end interface
    
    character :: clinoform_filename_P*1000,B*5, A*5, temp_name*1000, filepath*1000, endpath*1000
    integer :: SPACE, max, P, n, m, n_regrid, m_regrid, numsurf, i, j, k, n_regrid_2, m_regrid_2
    parameter (SPACE=6100)
    real :: nonvalue,xmin,xmax,ymin,ymax
    real, allocatable :: CLINOFORM_P(:,:,:)
    real, allocatable :: CLINOFORM_P2(:,:,:)
    real, allocatable :: CLINOFORM(:,:)
    real, allocatable :: passage(:,:)
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: NODE_CHANGE(:,:)
    real, allocatable :: NODE_CHANGE_2(:,:)
    
    allocate (passage(SPACE,5))
    allocate (TEMP(2,SPACE,SPACE)) 
    allocate (CLINOFORM_P(1,SPACE,SPACE))
    allocate (CLINOFORM_P2(1,SPACE,SPACE))
    allocate (CLINOFORM(SPACE,SPACE))
    allocate (NODE_CHANGE(SPACE,SPACE))
    allocate (NODE_CHANGE_2(SPACE,SPACE))
    
    call listing(SPACE,max)
    
    write(*,*)"Enter the file path for the heterogeneity surfaces"
    read(*,*) filepath 
    
    write(*,*)"Enter the file path destination for the heterogeneity surfaces"
    write(*,*)"This must be different to the file path for the input"
    write(*,*)"heterogeneity surfaces"
    read(*,*) endpath
    
do P=1,max
    call clinoform_filename(clinoform_filename_P, P, max, B, filepath)
    call clinoform_surfaces(TEMP,CLINOFORM_P,passage,clinoform_filename_P, n, m, P, nonvalue,xmin,xmax,ymin,ymax) 
    do i=1,n
      do j=1,m
          CLINOFORM(i,j)=CLINOFORM_P(1,i,j)
       end do
    end do
    !call more_nodes(CLINOFORM, n, m, n_regrid, m_regrid, NODE_CHANGE, Nonvalue, numsurf)
    !call more_nodes(NODE_CHANGE, n_regrid, m_regrid, n_regrid_2, m_regrid_2, NODE_CHANGE_2, Nonvalue, numsurf)
    call listing_surfaces(max,A,P)
    open(unit=4, file="listing_surfaces.txt")  
    close(unit=4, STATUS='DELETE')
    !do i=1,n_regrid_2
       !do j=1,m_regrid_2
         ! CLINOFORM_P2(1,i,j)=NODE_CHANGE_2(i,j)
       !end do
   ! end do
    k=1
    temp_name=trim(endpath)//"Heterogeneity_Surface"//trim(A)//".zmap+"
    call write_zmap_file(CLINOFORM_P,temp_name,n, m,k,xmin,xmax,ymin,ymax)
    write(*,*) "Completed Heterogeneity_Surface"//trim(A)//".zmap+"  
end do 

    open(unit=2, file="listing.txt")  
    close(unit=2, STATUS='DELETE')
    
    deallocate (passage)
    deallocate (TEMP) 
    deallocate (CLINOFORM_P)
    deallocate (CLINOFORM_P2)
    deallocate (CLINOFORM)
    deallocate (NODE_CHANGE)
    deallocate (NODE_CHANGE_2)
    
    
    write(*,*) 'Press return to exit program'
    read(*,*)
    
stop
end program heterogeneity_surfaces_ordered
!-------------------------------------------------------------------------------------------------------------------------------
subroutine array_allocation (temp_name,passage,TEMP,n,m, nonvalue,xmin,xmax,ymin,ymax)

    character :: temp_name*1000
    character :: dummy1*10, dummy2*10, dummy3*10, dummy4*10
    character :: dummy5*10, dummy6*10, dummy7*10, dummy8*10
    integer :: m, n, lines, rest, i, j, h, raw, mode, k, SPACE
    parameter (SPACE=6100)
    real :: numb, nonvalue, xmin, xmax, ymin, ymax, depth_min, depth_max
    real :: passagefive(5), passagefour(4), passagethree(3), passagetwo(2)
    real :: passageone
    real, allocatable :: passage(:,:)
    real, allocatable :: TEMP(:,:,:) 
    
   open(unit=1, file=temp_name)
	 read(1,*) dummy1
	 read(1,*) dummy2
	 read(1,*) dummy3
	 read(1,*) dummy4
	 read(1,*) dummy5
	 read(1,*) dummy6
	 read(1,*) numb, nonvalue
	 read(1,*) m, n, xmin, xmax, ymin, ymax
	 read(1,*) dummy7
     read(1,*) dummy8
     
     do i=1,SPACE
       do j=1,SPACE
          TEMP(1,i,j)=nonvalue
          TEMP(2,i,j)=nonvalue
       end do
    end do
      
     lines=m/5+1
     rest=m-(lines-1)*5
	if (rest==0) then;
    	lines=lines-1
    end if
    
	 do i=1,n
    	h=1
        do j=1,(lines-1)
           read(1,*) passagefive
           passage(h,1:5)=passagefive
           h=h+1
    	end do

		if (rest==1) then;
          read (1,*) passageone
          passage(h,1)=passageone
        else if (rest==2) then;
          read (1,*) passagetwo
          passage(h,1:2)=passagetwo
        else if (rest==3) then;
          read (1,*) passagethree
          passage(h,1:3)=passagethree
        else if (rest==4) then;
          read (1,*) passagefour
          passage(h,1:4)=passagefour
        else;
          read (1,*) passagefive
          passage(h,1:5)=passagefive
		end if
        
		do k=1,m
          raw=k/5
          mode=k-raw*5
          if (mode==0) then;
            mode=5
            raw=raw-1
          end if
          TEMP(1,i,k)=passage(raw+1,mode)          
         end do
     end do
    close(unit=1) 
              
return
end subroutine array_allocation
!----------------------------------------------------------------------------------------------------
subroutine listing(SPACE,max)

 integer :: i, max, SPACE

 open(unit=2, file="listing.txt")
  do i=1,SPACE
   write(2,*) i
  end do
 close(unit=2)
 
 write(*,*)"How many mud surfaces are there?"
 read(*,*) max

 return
end subroutine listing
!----------------------------------------------------------------------------------------------------
subroutine clinoform_filename(clinoform_filename_P, P, max, B, filepath)

    character :: clinoform_filename_P*1000,B*5, filepath*1000 
    integer :: i, P, max
          
    open (unit=2,file="listing.txt")
        do i=1,P
	       read(2,*) B
        end do
    close(unit=2)  
    clinoform_filename_P=trim(filepath)//"Heterogeneity_Surface"//trim(B)//".zmap+"
        
return
end subroutine clinoform_filename
!----------------------------------------------------------------------------------------------------
subroutine clinoform_surfaces (TEMP,CLINOFORM_P,passage,clinoform_filename_P, n, m, P, nonvalue,xmin,xmax,ymin,ymax)

   interface
     subroutine array_allocation (temp_name,passage, TEMP,n,m,nonvalue,xmin,xmax,ymin,ymax)
     character :: temp_name*1000
     integer :: n,m
     real :: nonvalue,xmin,xmax,ymin,ymax
     real, allocatable :: TEMP(:,:,:)
     real, allocatable :: passage(:,:)
     end subroutine
    end interface
    
   interface
     subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax, filepath)
     character :: temp_name*1000, filepath*1000
     integer :: n,m, k
     real :: xmin,xmax,ymin,ymax
     real, allocatable :: TEMP(:,:,:)
     end subroutine
    end interface

    character :: clinoform_filename_P*1000,temp_name*1000
    integer :: n, m, i, j, k, P, max
    real :: nonvalue,xmin,xmax,ymin,ymax
    real, allocatable :: passage(:,:)
    real, allocatable :: TEMP(:,:,:)
    real, allocatable :: CLINOFORM_P(:,:,:)
    
    temp_name=clinoform_filename_P
    call array_allocation (temp_name,passage,TEMP,n,m, nonvalue,xmin,xmax,ymin,ymax)
     do i=1,n
        do j=1,m
           CLINOFORM_P(1,i,j)=TEMP(1,i,j)
        end do
     end do
     !k=1
     !temp_name="CLINOFORM_test.zmap+"
     !call write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)
     !print*, "DUMMY"
         
return
end subroutine clinoform_surfaces
!----------------------------------------------------------------------------------------------------
subroutine write_zmap_file(TEMP,temp_name,n,m,k,xmin,xmax,ymin,ymax)

    character :: temp_name*1000
    character :: path*150
    integer :: n, m, i, k
    real :: xmin, xmax, ymin, ymax
    real, allocatable :: TEMP(:,:,:)
          
    
     open(unit=1,file=temp_name)
        write(1,*)'!'
        write(1,*)'!     ZIMS FILE NAME :  Clinoform_Surface.txt'
        write(1,*)'!     FORMATTED FILE CREATION DATE: xx/xx/xxxx'
        write(1,*)'!     FORMATTED FILE CREATION TIME: xx:xx'
        write(1,*)'!'
        write(1,1)'@Grid HEADER, GRID, 5'
1       format(A21)
        write(1,*)'	15,-9999.000,  ,	4 ,1.0000'
        write(1,"(I6,A,I5,A,1F12.4,A,1X,1F12.4,A,1F12.4,A,1F12.4)") m,',',n,',',xmin,',',xmax,',',ymin,',',ymax
        write(1,*)'0.0000     0.0000     0.0000'
	    write(1,2)'@'
2       format(A1)
        do i = 1, n 
            write(1,"(1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4)") TEMP(k,i,1:m)
         end do
    close(unit=1)    

return
end subroutine write_zmap_file
!-------------------------------------------------------------------------------------------------------------------------------
subroutine more_nodes (BIGbis, nbis, mbis, n_regridbis, m_regridbis, NODE_CHANGEbis,Nonvalue, numsurf)
 
 implicit none 
 
 integer :: nbis, mbis, cv, tv, SPACE, n_regridbis, m_regridbis, i, j, numsurf, choice, NODE
 parameter (SPACE=6100)
 parameter (cv = 2)
 parameter (tv = 2)
 parameter (NODE =6100)
 real, allocatable :: BIGbis(:,:)
 real, allocatable :: NODE_CHANGEbis(:,:)
 real :: Nonvalue
 
 
 
 n_regridbis = (2*nbis)-1
 m_regridbis = (2*mbis)-1
  
 do i = 1, n_regridbis
   do j = 1, m_regridbis
   
   if (MOD(i,2)==1.and.MOD(j,2)==1) then; 
        if(BIGbis(((i/2)+1),(j/2)+1)-1<Nonvalue) then;
        NODE_CHANGEbis(i,j)=Nonvalue
        else;
      NODE_CHANGEbis(i,j) = BIGbis(((i/2)+1),((j/2)+1))
      end if 
   else;
      
        if (MOD((i+j),2)==0) then; 
            if(BIGbis(((i/2)+1),(j/2)+1)-1<Nonvalue.or.BIGbis((i/2)+1,(j/2))-1<Nonvalue.&
            &or.BIGbis((i/2)+1,(j/2))-1<Nonvalue.or.BIGbis((i/2),(j/2)+1)-1<Nonvalue) then;
            NODE_CHANGEbis(i,j)=Nonvalue
            else;
            NODE_CHANGEbis(i,j) = (BIGbis(((i/2)+1),(j/2)+1)+BIGbis((i/2)+1,(j/2))+BIGbis((i/2)+1,(j/2))+BIGbis((i/2),(j/2)+1))/4 
            end if       
        else;
           
           if (MOD(i,2)>0) then;
                if(BIGbis(((i/2)+1),(j/2)+1)-1<Nonvalue.or.BIGbis((i/2)+1,(j/2))-1<Nonvalue) then;
                NODE_CHANGEbis(i,j)=Nonvalue
                else; 
                NODE_CHANGEbis(i,j)=(BIGbis((i/2)+1,(j/2)+1)+BIGbis((i/2)+1,j/2))/2
                end if
           else;
                if(BIGbis((i/2)+1,(j/2)+1)-1<Nonvalue.or.BIGbis((i/2),(j/2)+1)-1<Nonvalue) then;
                NODE_CHANGEbis(i,j)=Nonvalue
                else;
                NODE_CHANGEbis(i,j)=(BIGbis((i/2)+1,(j/2)+1)+BIGbis((i/2),(j/2)+1))/2
                end if
           end if
        
        end if
        
    end if
    end do
  end do
  
  
 return
 
 end subroutine more_nodes
 !-------------------------------------------------------------------------------------------------------------------------------
 subroutine listing_surfaces(max,A,P)
 
 character :: A*5
 integer :: g,max,P
 
 g=max-(P-1)
 
 open(unit=4, file="listing_surfaces.txt")
   write(4,*) g
   !print*, g
 close(unit=4)
 
 open(unit=4, file="listing_surfaces.txt")
   read(4,*) A
 close(unit=4)

 return
end subroutine listing_surfaces
!----------------------------------------------------------------------------------------------------